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Abstract — There  has  been  an  increasing  interest  in  on-orbit  au¬ 
tonomous  servicing  and  repair  of  satellites  as  well  as  controlled 
active  debris  removal  (ADR)  in  the  space  industry  recently.  One 
of  the  most  challenging  tasks  for  servicing/repair  as  well  as  for 
ADR  is  the  rendezvous  and  docking  with  a  non-cooperative 
tumbling  resident  space  object  (RSO).  This  paper  presents  a 
propellant  optimal  maneuver  profile  for  a  servicing  spacecraft 
to  perform  proximity  operations  and  eventually  dock  with  a 
non-cooperative  target.  The  strategy  is  to  find  an  optimal 
trajectory  which  will  guide  the  servicing  spacecraft  to  approach 
the  tumbling  satellite  such  that  the  two  vehicles  will  eventually 
have  no  relative  motion.  Therefore,  a  subsequent  docking  or 
capture  operation  can  be  safely  performed.  The  research  de¬ 
scribed  here  elaborates  on  the  previous  work  that  studied  the 
minimum-control-effort  for  a  3-D  docking  to  a  tumbling  object 
considering  a  full  six-degree-of-freedom  model  of  both  chaser 
and  target.  The  current  work  expands  the  scope  by  adding  new 
set  of  linearized  equations  of  motion  that  capture  the  effect  of 
the  J2  geopotential  disturbance  force. 

Typically,  Hill’s  linearized  equation  of  relative  motion  have  been 
used  for  this  analysis,  but  they  fail  to  capture  the  effect  of  J2 
disturbance  force  on  the  chaser  satellite.  Firstly,  the  effects  of 
the  J2  disturbance  force  is  added  to  the  linearized  equations  of 
motion  by  the  addition  of  the  J2  terms.  Secondly,  minimum- 
control-effort  optimality  condition  is  examined  and  propellant 
optimal  trajectories  for  a  relative  motion  problem  are  then 
numerically  solved,  by  using  a  direct  collocation  method  based 
on  the  Gauss  pseudospectral  approach.  The  simulation  results 
shows  the  effect  of  errors  caused  by  the  oblateness  of  the  earth 
(as  described  by  the  J2  potential)  on  the  described  relative 
motion  problem.  Furthermore,  effect  of  J2  disturbance  on  the 
optimal  trajectory  is  discussed  for  the  minimum  propellant- 
consumption  optimality  condition. 
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1.  Introduction 

There  has  been  an  increasing  interest  in  satellite  on-orbit 
autonomous  servicing  in  the  space  industry  recently.  JAXA 


recently  completed  a  technology  demonstration  mission  ETS- 
7  [13].  NASA  did  an  autonomous  rendezvous  mission 
through  the  Demonstration  for  Autonomous  Rendezvous 
Technology  (DART)  in  2005  [8],  where  the  mission  was  not 
completed  due  to  more  than  expected  propellant  usage  during 
rendezvous  maneuvering.  Air  Force  Research  Laboratory’s 
(AFRL)  Experimental  Satellite  Systems-10  and  11  (XSS-10 
and  XSS-1 1)  have  been  developed  in  order  to  show  the  ability 
for  a  small  sat  to  autonomously  plan  and  rendezvous  with 
a  passive  and  cooperative  Resident  Space  Object  (RSO)  in 
Low  Earth  Orbit  (LEO)  [1].  In  addition.  Defense  Advance 
Research  Projects  Agency’s  (DARPA)  Orbital  Express  (OE) 
Advance  Technology  Demonstration  Program  validated  the 
technology  and  techniques  for  on-orbit  refueling  and  configu¬ 
ration  of  two  satellites  [3]  [4].  The  existence  of  this  programs 
demonstrates  that  there  is  a  need  for  a  robust  and  effective 
autonomous  close  proximity  control  algorithms  for  multiple 
spacecraft.  The  use  of  micro-satellites  to  inspect,  service, 
repair,  deorbit  and  refuel  larger  spacecrafts  is  a  long-term 
goal. 

In  order  to  perform  on-orbit  service,  the  servicing  spacecraft 
has  to  first  rendezvous  and  dock  with  the  satellite  to  be 
serviced  in  orbit.  The  tumbling  of  an  uncontrolled  resident 
space  object  (RSO)  or  the  rotational  and  transnational  motion 
of  an  non-cooperating  RSO  present  challenges  to  the  docking 
operations. 

From  a  theoretical  standpoint,  the  present  paper  elaborates 
on  the  previous  work  by  Boyarko  et  al.  [5][6][7],  who 
studied  the  minimum-time  and  minimum-control  problem 
for  a  rendezvous  of  a  chaser  satellite  to  a  tumbling  object 
considered  a  full  six-degree-of-freedom  model  of  both  chaser 
and  target.  The  current  work  expands  the  scope  by  adding  the 
effect  of  J2  geopotential  disturbance  force  and  investigates  its 
effects  on  the  proximity  operations. 

The  majority  of  previous  research  pertaining  to  operations 
that  focused  on  rendezvous  with  an  uncontrolled  RSO,  the 
relative  motion  dynamics  pertinent  to  proximity  space  opera¬ 
tion  only  considers  Hill’s  equations,  also  known  as  Clohessy- 
Wiltshire  equations.  The  reason  being  Hill’s  equations  are 
very  simple  to  implement  and  have  been  successfully  used 
to  describe  the  relative  motion  of  two  satellite  during  ren¬ 
dezvous  maneuvers.  Although,  the  assumption  made  by 
Hill’s  equations  that  is  cited  repeatedly  as  the  main  source 
of  error  is  that  the  Earth  is  perfectly  spherical.  Because  the 
Earth  is  not  perfectly  spherical,  but  rather  an  oblate  spheroid 
and  the  dominant  term  of  the  series  expansion  of  this  effect  is 
the  J-2  term,  modeling  errors  are  introduced  by  the  noncentral 
forces.  To  reduce  the  effect  of  the  modeling  errors  that  are 
present  in  Hill’s  equations,  much  research  has  been  done  with 
varying  success  to  incorporate  the  effects  of  the  J2  potential. 
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This  papers  implements  the  formulation  of  Schweighart  and 


Sedwick  [11]  [12]  which  describes  the  linear  combination  of 
Hill’s  equation  with  J2  effect.  The  equations  derived  by 
Schweighart  and  Sedwick  [11]  [12]  provides  insight  into  the 
relative  motion  of  satellites  under  the  influence  of  the  J2  po¬ 
tential.  The  period  mismatch  between  the  in-plane  and  cross¬ 
track  motion  produces  an  effect  called  tumbling  in  which 
the  spacecraft  formation  appears  to  tumble  around  the  orbital 
angular  momentum  vector.  Additionally,  small  differences  in 
the  orbital  inclination  of  the  satellites  in  the  formation  cause 
differential  drift  in  the  longitude  of  the  ascending  node.  This 
differential  drift  causes  the  spacecraft  formation  to  break  up 
over  time.  The  paper  is  organized  as  follows.  Section  2 
presents  the  equations  of  relative  motion  that  incorporate  the 
effect  of  the  J2  perturbation.  The  analysis  of  the  relative 
motion  with  Hill’s  equations  and  the  derived  time-averaged 
J2  model  with  corrections  along  with  the  full  J2  effect  is 
presented  in  Section  3.  The  set  up  of  the  optimization 
problem  is  presented  in  Section  4  and  the  results  are  presented 
and  analyzed  in  Section  5. 


2.  Relative  Motion  Models 

Figure  [1]  illustrates  the  curvilinear  x  (radial)  y  (along-track) 
z  (cross-track)  axis  system  used  for  the  problem.  The  x  vector 
points  in  the  radial  direction,  the  z  vector  is  perpendicular  to 
the  orbital  plane  and  points  in  the  direction  of  the  angular 
momentum  vector.  Finally,  they  vector  completes  movement. 
In  the  x  -  y  -  z  coordinate  system,  the  specification  is  made 
that  it  is  a  curvilinear  coordinate  system.  The  x  vector 
remains  unchanged,  however  the  y  and  the  z  vector  ’curves’ 
around  the  orbit. 

Unlike  the  Local  Vertical,  Local  Horizontal  (LVLH)  body 
fixed  frame,  implemented  by  the  Hill’s  equations,  the  only 
difference  comes  from  the  fact  that  the  LVLH  frame  is  not 
defined  as  a  curvilinear  system  but  as  a  rectangular. 


n 


Figure  1.  The  x-y-z  Curvilinear  Body-Fixed  Coordinate 
System. 

Figure  [2]  illustrates  the  r  -  6  -  i  spherical  coordinate  system 
used  to  describe  the  gravitational  acceleration  and  the  J2 
potential  due  to  the  spherical  earth  .  The  r  points  in  the  radial 
direction,  i  is  the  azimuthal  angle  measured  around  the  line  of 
nodes,  and  6  is  the  co-latitude  measured  from  the  ascending 
node  which  acts  as  the  pole  of  the  spherical  system. 
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Figure  2.  The  r-6-i  Spherical  Body-Fixed  Coordinate  Sys¬ 
tem. 


Conversion  between  both  these  body  fixed  coordinate  sys¬ 
tems  is  fairly  straightforward.  The  radial  vector  in  each 
coordinate  system  completely  coincides.  The  y  vector,  and 

6  coincide;  and  the  z  vector  and  the  i  vector  coincide.  In  this 
way,  the  curviliniar  coordinate  system  defined  is  very  much 
like  a  spherical  coordinate  system. 

A  summary  of  the  mathematical  model  development  is  per¬ 
formed  and  described  below  in  order  to  highlight  the  set 
of  equations  of  motion  that  are  considered  for  the  relative 
motion  problem  between  a  controlled  satellite  (Chaser)  w.r.t 
to  the  assumed  reference  orbit  (Non-Cooperative  Target). 

The  derivation  begins  with  the  analytic  equation  of  motion  for 
an  orbiting  satellite  under  the  influence  of  the  J2  potential. 

f  =  g(r)+J2{r)  (1) 

where,  g(r)  is  the  gravitational  acceleration  due  to  spherical 
Earth, 

g(r)  =  -( n/r2)r  (2) 

J 2  (r* )  is  the  acceleration  due  to  J2  potential  [9], 


/2M  =  -(3/2)(J2[J,Rl/r4)l(l  -  3sin2  isin2  6)x 
+(2  sin2  i  sin  6  cos  0)y  +  (2  sin  i  cos  i  sin  6)z] 


r  is  the  position  vector  of  the  satellite,  and  x-y-z  is  the 
coordinate  system  described  in  Figure  [1]. 

The  models  are  developed  by  considering  the  relative  motion 
between  one  satellite  and  a  reference  orbit.  For  now,  the 
derivation  only  requires  that  the  reference  orbit  is  constant 
radius  (circular)  and  the  equation  of  motion  for  the  reference 
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(1  -  3  sin2  i  sin2  9) 
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sin2  i  sin  29 
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orbit  is  given  by  the  Newton's  law  of  motion  in  a  central  force 

field:  x  =  r  -  rlei  -  2lj  x  x  -  u  x  x  -  u  x  (u  x  x)  (9) 


»W  =  gOref)  (4) 

Furthermore,  linearizing  the  gravitational  terms  in  Equation 
(1)  with  respect  to  the  defined  reference  orbit  results  in 

r  =  gOref)  +  g(fref)  x+J  2  (rref)  +  VJ2  (rre f)  •  X  (5) 


For  a  circular  reference  orbit  the  rotation  rate  of  the  coordi¬ 
nate  system  is  constant  given  by 


w  =  nz, 


(10) 


where,  the  relative  position  of  the  satellite  with  respect  to  the 
reference  orbit  is  given  by  x  as  represented  in  Figure  [3] 


Substituting  Equation  (5)  into  Equation  (9)  and  rearranging 
the  terms  gives  us: 


x  =r  -  rref 


(6) 


x  +  2uxx+wxx+ux  (u  x  x)  =  g(rref) 
+g(fref)  •  X  +  J2  (r ref)  +  V/2  (fref)  ’  X  -  rIe f 
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The  equation  above  can  be  further  simplified  by  canceling 
out  g(rref)  with  rlef  as  they  are  equal  as  seen  in  equation 
(4).  In  the  following  sections  the  above  equation  along 
with  Hill’s  Eqn.  are  presented  as  a  comparison  to  the  time 
averaged  J2  equations,  with  reference  orbit  period  and  nodal 
drift  corrections.  These  time  averaged  J2  equations  with 
corrections  are  derived  in  the  following  subsections  below. 

Time  Averaged  V  J2 

Equation  (11)  is  a  linearized  differential  equation  of  motion 
with  time  varying  coefficients.  The  ODE  with  time  varying 
coefficient  in  Equation  (11)  is  cast  into  an  ODE  with  time 
varying  coefficient  by  taking  the  time  average  of  the  VJ2 
term 


1  [2\j2(r)d9=  4 

As  0 

0  ‘ 

Figure  3.  The  Reference  Orbit 
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When  a  spherical  coordinate  system  (r  -  9  -  i )  with  the  pole 

aligned  with  the  ascending  node  is  used,  the  gradient  of  the 
g(r)  gravitational  acceleration  can  be  calculated.  The  result 

where. 

is: 

s  =  (3  J2f?2/8r2)(l  +  3  cos  2i) 

03) 

Vg(r) 


2  {n/r3)  0  0 

0  — (/r/r3)  0 

0  0  -(ft/r3) 


Adding  the  time-averaged  term  and  removing  g(rref)  with  rref 
Equation  (11)  becomes, 


The  J2  disturbance  in  Equation  (3)  is  given  in  i  -  y  - 
z  coordinates.  However,  the  equation  can  be  transformed 

directly  to  the  r  -  6  -  i  coordinate  system  without  any  loss 
of  generality.  Taking  the  gradient  gives  Equation  (8) 


x  +  2uixx+u)xx+lox  (w  x  x)  =  g{rref)  ■  x 

1  f2K  04) 

+J2(rref)  +  —  J  V/2(rref)  ■  xd9 


Because  the  reference  frame  is  rotating,  rotational  terms  are 
added  when  calculating  the  relative  acceleration  and  veloci¬ 
ties  of  the  satellite  w.r.t.  the  reference  orbit. 
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By  adding  the  time  average  of  the  gradient  of  the  J2  distur¬ 
bance,  the  periodic  component  of  the  gradient  is  lost.  As  the 
values  of  g(rref)  and  /2(rref)  and  the  gradients  of  these  terms 


are  evaluated  at  the  reference  orbit,  x  must  remain  relatively 
small  for  the  the  first  order  term  of  the  Taylor  series  expansion 
that  is  used  in  the  derivation  of  the  equations  of  motion  to 
apply  and  therefore,  the  mean  radii  and  orbital  rates  of  the 
perturbed  satellite  and  reference  satellite  must  be  similar. 

Effect  of  the  J-2  perturbation  also  causes  a  drift  in  the  longi¬ 
tude  of  ascending  node  which  varies  with  inclination.  The  J2 
perturbed  satellite  will  experience  this  drift,  but  the  reference 
orbit  will  not,  thus  the  two  orbits  will  gradually  separate 
unless  some  correction  to  the  orbital  rate  and  nodal  drift  of 
the  reference  satellite  is  included  in  the  model. 

The  fact  that  the  VJ2  term  has  been  time  averaged  over  a 
period  in  Equation  (14)  has  the  effect  to  cause  errors  in  cross¬ 
track  motion.  It  is  seen  that  cross-track  motion  is  due  solely 
to  that  the  satellite’s  orbit  and  the  associated  reference  are 
not  coplaner.  It  is  a  periodic  motion  that  is  equal  to  zero 
when  the  two  orbital  planes  intersect  and  is  at  a  maximum 
90°  away  from  the  intersection  of  the  planes.  As  this  paper 
deals  with  docking  application  focusing  where  the  satellite  is 
in  close  vicinity  with  the  reference  orbit,  both  having  similar 
inclinations,  effects  to  the  cross-track  errors  are  not  accounted 
in  this  paper.  Although  addition  and  studying  the  effect  of 
these  cross  track  errors  will  be  a  topic  explored  in  future. 

Correcting  the  orbital  period  of  the  reference  orbit 

Under  the  influence  of  the  J2  disturbance,  the  perturbed 
satellite  will  have  a  different  orbital  period  than  when  unper¬ 
turbed.  Because  of  this  discrepancy,  the  perturbed  satellite 
drifts  from  the  unperturbed  reference  orbit  and  eventually  the 
linearized  equations  break  down.  To  fix  this  problem,  the 
period  of  the  reference  orbit  must  be  adjusted  to  match  the 
period  of  the  satellites  in  the  cluster. 

The  change  in  period  due  to  the  J2  disturbance  as  given  by 
Schweighart  and  Sedwick  can  be  found  from  the  average  J2 
disturbance  (not  to  be  confused  with  the  time  average  of  the 
gradient  of  the  J2  term  taken  above).  The  equation  of  motion 
of  the  reference  orbit,  as  seen  in  Equation  (4),  now  becomes, 

1  [2n 

Uef  =  g(rKf)  +  —  /  J-2(rK{)d6  (15) 

47r  J q 


As  the  rate  of  the  reference  orbit  has  changed,  so  has  the 
average  angular  speed  of  the  reference  satellite,  and  the 
coordinate  system  which  is  based  there.  The  new  angular 
velocity  can  be  found  by  equating  the  accelerations  to  give 

u  =  ncz,  c  =  Vl+~s  (18) 

where  ”s”  in  Equation  (15)  and  Equation  (17)  is  given  by, 

a  =  3  (1  +  cos  2 i)  (19) 

8rz 

Correcting  the  Reference  Orbit  for  Nodal  Drift 

Although  the  preceding  equations  of  motion  are  a  vast  im¬ 
provement  over  Hill’s  equations  when  incorporating  the  J2 
disturbance,  more  can  be  done.  Even  though  the  orbital 
period  of  the  reference  orbit  has  been  adjusted  to  match  the 
period  of  the  satellite  under  the  influence  of  the  J2  poten¬ 
tial,  they  still  drift  apart  due  to  separation  of  the  longitude 
of  the  ascending  node.  Schweighart  and  Sedwick  worked 
out  analytically  a  separation  distance  between  the  perturbed 
satellite  and  the  reference  orbit  as  rref  AO  sin  i  over  one 
period.  Similarly,  after  two  orbital  periods,  the  satellite  and 
the  reference  orbit  will  be  separated  by  rref2AO  sin  i,  and  this 
process  will  continue  to  cause  them  to  drift  farther  and  farther 
apart.  Although  the  equations  of  motion  in  their  current  form 
do  capture  the  bulk  of  this  motion,  the  increasing  separation 
causes  the  linearization  to  break  down  after  several  periods. 
Because  of  this,  the  reference  orbit  must  again  be  redesigned 
so  that  it  tracks  this  secular  motion. 

This  separation  in  longitude  of  ascending  node  as  given 
by  Schweighart  and  Sedwick,  shows  that  only  the  normal 
component  of  the  J2  acceleration  is  responsible  for  the  drift. 
Applying  the  normal  component  of  the  J2  disturbance  accel¬ 
eration  to  the  updated  reference  orbit  [Equation  (15)]  gives 
us: 


1  f2’7 

Aef  =  g(fm f)  +  ^  J  J2{r ref) dO  +  [J2(r ref)  '  z]z  (20) 


where  the  time  averaged  J2  acceleration  is  given  by  Eq.  (16) 
below. 


1 

2-7T 


/•2tt 

—n2rs 

/  Mr)d0  = 

0 

1 0 

0 

(16) 


This  results  in  a  non-zero  value  for  acceleration  in  the  radial 
x  direction  only.  This  can  be  visualized  as  an  additional  force 
acting  to  increase  the  Keplerian  gravity  term.  Therefore,  if  a 
satellite  is  to  remain  in  a  circular  orbit  its  orbital  rate  must  be 
increased  above  that  for  Keplerian  dynamics. 

Substituting  Equation  (15)  in  Equation  (11),  we  get 


X  +  ZutXX+UXX+UX  (u  XX)  =  g{r ref)  X  +  J2{rIef) 


1  f27T 

+  —  J  Vy 2  (f ref)  d9  ■  X 


1 

27T 


Time-Averaged  J2  EQM  with  Corrections 

The  final  equation  of  motion  which  captures  the  effect  of  J2 
disturbance  alonga  with  reference  orbit  period  and  nodal  drift 
corrections  is  described  below.  As  we  added  a  component 
of  the  J2  potential  to  the  reference  orbit  for  the  correction 
of  nodal  drift,  that  component  must  be  subtracted  from  the 
equations  of  motion  of  the  satellite  relative  to  the  reference 
orbit  [Equation  (17)].  The  resulting  equations  of  motion  in 
vector  for  are: 


X  +  2UXX+U1XX+U1X  (wxi)=  g(r  ref)  ■  x 
1  f27T 

A  J  2  (j" ref )  A  ~  /  S7J2(r ref)  dO  ■  X 

Z7T  J 0 

1  f271 

~7T  J2(rref)d0  -  [J2{rie f)  •  z]z 

Jo 


J  2  0*ref  )d0 
(17) 
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The  above  equations  can  be  expanded  by  substituting  the 
appropriate  terms  results  in  a  set  of  differential  equations  of 
form. 


x  —  2  {nc)y  —  (5c2  —  2  )n2x  =  —3n2J2(R2/rrei) 
x  [|  —  3  sin2  rref  sin 2(kt) /2  —  (1  +  cos  2zref )  / 8] 


(22) 


y  +  2(nc)x  =  —3n2J2(R2/2rref)  sin2  iref  sin(2fcf) 
z  +  (3c2  —  2)n2z  =  0 


where  k  is  defined  as. 


k  =  nc  + 


3^jU2R2e 

2r7/2 


cos2  i 


(23) 


The  above  differential  equations  are  used  for  solving  the 
optimal  control  problem  described  in  this  paper.  Prior  to 
setting  up  the  simulation,  the  section  below  goes  over  the 
effect  of  J2,  with  and  without  corrections  as  compared  to 
Hill’s  equation,  to  highlight  the  motivation  of  this  study  and 
understand  the  effect  of  the  errors  on  the  relative  position  of 
the  perturbed  satellite. 
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3.  Effects  of  J2  disturbance 

This  section  highlights  the  effect  of  .J2  perturbations  on  the 
relative  motion  of  the  perturbed  satellite  with  respect  to  the 
reference  orbit.  Moreover,  error  related  to  reference  orbital 
period  and  nodal  drift  are  highlighted.  Finally,  the  effect  of 
disturbance  of  J2  potential  is  studied  in  comparison  to  Hill’s 
equation  of  motion. 

Numerical  simulations  are  performed  using  fixed-step  ODE 5 
integrator  to  plot  relative  motion  trajectories.  The  orbital 
parameters  chosen  for  the  reference  orbit  where  the  target 
is  placed  are:  ire f  =  70°,  the  apogee  and  perigee  altitude 
is  given  by  ra  =  919-km  and  rp  =  902-km  respectively. 
The  rationale  for  choosing  the  above  orbital  parameters  is 
discussed  in  Sec.  4  where  the  optimization  problem  for  this 
study  is  formulated. 

The  following  Figure  [4]  demonstrates  the  effect  of  the  J2 
geopotential  disturbances  at  70°  inclination.  The  purpose  to 
acquire  the  J2  effect  on  the  chaser  from  both  the  equatorial 
and  polar  regions,  a  logical  choice  of  starting  from  [-1,  0, 
0]  km  is  employed  as  the  initial  condition  for  the  relative 
position,  thus  in  order  to  attempt  to  only  simulate  the  J2  effect 
on  the  pure  in-track  motion.  The  chaser  is  assumed  to  be 
at  rest  initially.  All  the  simulations  below  for  this  section 
are  propagated  for  the  time  span  of  two  orbital  period  of  the 
reference  orbit  which  comes  to  3.44  hrs. 

It  is  seen  in  Figure  [4]  the  considerable  relative  position 
error  created  by  the  J2  disturbance  and  the  importance  for  its 
implementation  while  studying  the  optimal  trajectory  design 
for  rendezvous  and  proximity  operations.  Additionally,  it  is 
noted  that  the  J2  model  without  corrections  corresponding  to 
Equation  (11)  adds  sizable  cross-track  errors  as  a  result  of 
orbital  period  mismatch  between  the  target  and  the  chaser 
orbits  and,  the  nodal  drift  caused  by  the  separation  of  the 
longitude  of  the  ascending  node. 


error  in  the  in-track  motion  (x)  tops  out  at  ~  15  km  and  error 
propagation  in  time  looks  periodic.  Similarly,  for  the  radial 
motion  (y)  the  error  due  to  the  J2  seems  to  be  increasing 
sharply  with  time.  Finally,  for  the  cross-track  motion  (z) 
even  though  the  error  is  not  significant  compared  to  in-track 
and  radial  motion,  the  profile  shows  an  increasing  deviation 
with  time.  With  periodic  errors  in  the  in-tack  motion  and  the 
increasing  errors  in  both  the  radial  and  cross-track  direction, 
a  significant  effect  of  .J2  perturbation  on  the  chaser  is  seen. 


Time  [s] 


Figure  5.  Relative  position  errors  between  Hill’s  Eqns.  and 
Time-averaged  J2  with  corrections  Eqns.  (iref  =  70°). 


The  relative  position  errors  between  the  commonly  used 
Hill’s  equations  and  the  time-averaged  J2  with  corrections 
model  is  shown  in  Figure  [5],  As  seen  in  the  figure,  the 
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Figure  [6]  presents  similar  relative  position  errors  between 
the  time-averaged  J2  model  with  corrections  and,  the  J2 
equations  without  any  corrections  [Equation  (11)]  to  high¬ 
light  the  effect  caused  by  the  orbital  period  mismatch  and 
the  separation  in  the  ascending  node  causing  the  nodal  drift. 
As  expected,  the  position  error  in  all  three  directions  is 
significant.  As  a  consequence,  the  implementation  of  the 
corrections  is  crucial  for  studying  the  optimal  trajectory  de¬ 
sign  for  proximity  applications.  Moreover  when  compared  to 
Figure  [5],  it  is  seen  that  the  implementation  of  the  correc¬ 
tions  reduces  error  in  the  z  direction  which  causes  the  motion 
of  the  chaser  to  create  heavy  cross-track  deviations  from  the 
nominal  as  seen  in  Figure  [4] . 


Figure  6.  Relative  position  errors  between  Time-averaged 
J2  with  corrections  Eqns.  and  J2  without  corrections  Eqns. 
[Eqn.  (ll)](w  =  70°). 
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Figure  7.  Relative  motion  trajectories  (iref  =  0°). 


An  interesting  study  would  be  to  analyze  the  same  three 
models  at  0°  inclination  to  validate  the  physics  of  the  derived 
equations.  Figure  [7]  shows  the  relative  motion  trajectories 
under  the  influence  of  the  J2  disturbance  force.  Is  is  seen 
that  the  relative  motion  of  the  chaser  under  the  time-averaged 
J2  with  corrections  model  resembles  with  that  of  Hill’s  equa¬ 
tions,  as  logically  there  would  not  be  any  effect  of  J2  across 


the  equator.  Although  a  significant  error  is  seen  in  the  relative 
position  given  by  the  J2  equations  without  corrections. 

Figure  [8]  and  Figure  [9]  illustrates  the  errors  in  the  relative 
position  of  the  spacecrafts  between  J2  model  with  correc¬ 
tions  with  Hill’s  Equations  and  J2  model  without  corrections 
[Equation  (11)]  respectively. 
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Figure  8.  Relative  position  errors  between  Hill’s  Eqns.  and 
Time-averaged  J2  with  corrections  Eqns.  (zref  =  0°). 


Figure  9.  Relative  position  errors  between  Time-averaged 
J2  corrections  Eqns.  and  J2  without  corrections  Eqns.  [Eqn. 

(11)]  (iref  =  0°). 
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Figure  [8]  shows  that  the  position  errors  between  J2  cor¬ 
rections  model  with  Hill’s  equation  are  very  small.  The 
insignificant  error  noticed  corresponds  to  the  curvilinear  co¬ 
ordinated  system  used  for  the  derivation  of  the  time-averaged 
J2  equations.  Validating  the  time-averaged  J2  corrections 
model  once  again  at  0°  inclination  proves  that  the  model  used 
for  the  optimization  problem  formulation  in  the  following 
section  is  a  generalized  case  of  Hill’s  Equation  with  J2 
perturbation  incorporated  in  it. 

Figure  [9]  further  shows  the  inability  of  the  full  J2  model 
without  correction  [Equation  (11)]  to  fully  capture  the  orbital 
period  mismatch  and  nodal  drifts  errors.  As  expected  the 
position  errors  in  all  three  directions  are  significant. 


Using  the  controls  the  two  spacecrafts  are  brought  together 
from  some  initial  condition,  given  by  6  initial  value  of  states 
Xi(to),  i  =  1, 6  to  docking  conditions  described  by  a 
final  state,  which  for  this  problem  is  assumed  to  be  a  null 
i.e.  Xf  =  O. 

The  goal  of  the  optimization  problem  is  to  minimize  propel¬ 
lant  consumption  during  a  maneuver  of  finite  time  of  duration 
tf.  The  objective  function  is  then: 

f*f 

J{ u)  =  /  (|  Wi(f)  |  +  |  U2(t)  I  +  I  u3(t)  I )dt  (27) 

Jo 


From  the  study  conducted  in  this  section  we  see  the  im¬ 
portance  of  implementation  of  J2  in  the  linear  equations  of 
relative  motion  for  the  study  of  the  optimization  operations. 


4.  Optimization  Problem  Formulation 

This  section  develops  a  model  of  target-chaser  rendezvous. 
The  curvilinear  coordinate  system  described  in  Section  2  is 
used  for  the  model  formation.  We  start  from  the  arbitrary 
relative  position  and  would  like  to  bring  the  two  spacecraft 
together  for  docking. 

Using  the  notations  described  in  the  beginning  of  this  paper 
based  on  the  above  model,  the  dynamics  of  the  two  system 
can  now  be  described  as  follow.  The  translational  kinematics 
and  dynamics  of  a  chaser  spacecraft  in  the  orbit  frame  cen¬ 
tered  at  the  target  vehicle  are  given  by  Equation  (21), 


Optimization  Problem  Setup 

A  sample  maneuvering  scenario  is  considered  with  the  chaser 
center  of  mass  starting  at  a  distance  of  1  km  behind,  1  km 
below,  and  1  km  sideways  from  the  target  center  of  mass. 
Hence  the  initial  values  of  the  state  for  computer  simulations 
discussed  in  the  following  sections  is: 

X0  =  [-1000,  -1000,  -1000, 0, 0, 0] 

The  problem  is  further  simplified  by  neglecting  the  geometri¬ 
cal  structure  of  the  both  chaser  and  target.  For  the  time  being 
the  final  state  is  at  the  CoM  of  the  target  and  with  null  relative 
velocity  vector. 

X/  =  [0,0, 0,0, 0,0] 


x  =  2  (nc)y  +  (5c2  —  2  )n2x  —  3n2  J2(i?2/rref) 
x  [g  —  3  sin2  rref  sin2(fci) /2  —  (1  +  cos  2iref) /8]  +  fx 


y  =  —2 (nc)x  —  3n2  J2(i?2/2rref)  sin2  iref  sin(2fcf)  +  fy 


z  =  —(3c2  —  2  )n2z  +  fz 


(24) 


where  fx,  fy  and  fz  are  applied  forces  (controls)  expressed 
in  the  Hill’s  frame. 

Equation  (24)  hence  defines  a  6-state  system  of  differential 
equations  of  relative  motion  governing  the  docking  problem 
dynamics.  Combined  into  the  state  vector  X  these  states  are, 

X  =  [x,y,z,x,y,z\  (25) 


The  governing  dynamics  assumes  three  normalized  controls: 


u  = 


fx 


fy 

ymax 


fz 


(26) 


For  simplicity,  it  is  assumed  that  fimax  =  lm/s2  for  i  = 
x ,  y,  z.  Once  again,  these  controls  are  three  normalized 
components  of  a  translational  force  acting  on  a  chaser  /, 
( i  =  x,  y,  z ),  expressed  in  Hill’s  coordinate  frame.  All  three 
controls  are  bounded:  umjn  <  u  <  umax.  The  details  of 
the  thrust  bounds  are  discussed  in  the  subsection  along  with 
chaser  mass  and  its  thrust  parameters. 


Baseline  parameters  are  selected  from  a  proposed  controlled 
active  debris  removal  (ADR)  mission  Curimba  by  Udrea  and 
Nayak  [14],  The  authors  describes  Agena-D  RBs  as  one  of 
the  potential  target  for  the  mission.  Hence,  the  orbital  param¬ 
eters  for  this  upper  stage  debri  corresponds  to  the  reference 
orbit  chosen  for  this  problem.  As  most  of  the  Agena-D  RBs 
are  populated  at  approximately  70°  the  reference  inclination 
for  this  problem  is  assumed  70°.  Additionally,  the  apogee 
and  the  perigee  altitude  selected  for  the  reference  orbit  where 
our  target  spacecraft  is  placed  resembles  that  of  the  Agena-D 
RB  which  corresponds  to  ra  =  919  km  and  rp  =  902  km 
respectively.  Furthermore,  a  time  constraint  is  added  to  the 
simulation  by  fixing  the  final  time  to  one-orbital  period  of 
the  reference  orbit,  as  for  the  ADR  mission  it  requires  the 
chaser  to  rendezvous  with  the  target  with  the  shortest  amount 
of  time. 

The  mission  proposes  a  preliminary  design  of  the  chaser  by 
considering  a  12U  CubeSat.  This  CubeSat  is  a  three-axis 
stabilized  spacecraft  of  240x240x360  mm  and  a  mass  of 
24  kg.  It  is  known  that  the  CubeSat  has  a  miniaturized 
system  consisting  of  16  solenoid-fed  thrusters  placed  such 
that  they  generate  positive  and  negative  torques  about  all  three 
body  axes  with  dual  redundancy.  The  thruster  arrangement 
also  generates  positive  and  negative  forces  about  two  of  the 
three  body  axes.  For  the  purpose  of  this  work  it  has  been 
assumed  that  the  force  along  the  third  body  axis  is  available 
when  required.  Future  research  will  include  constraints  on 
the  generation  of  force  about  the  third  body  axis  to  take  into 
account  the  finite  slew  time. 

It  is  known  from  the  preliminary  design  provided  by  the 
authors  that  the  propellant  used  will  be  1,1-  difluorethane 
which  can  generate  about  70  s  of  specific  impulse.  As  for  the 
first  mission  design  iteration  they  assumed  that  each  thruster 
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produces  50  mN  of  thrust  force  [14].  Knowing  now  the  thrust 
force  for  individual  thrusters  (50  mN)  and  the  chaser  mass 
(24  kg),  control  bounds  can  be  calculated  for  the  optimal 
problem  for  the  Equation  (26).  By  these  means  the  control 
bounds  thus  will  be  calculated  by  ±FT/mchaser-  where  Ft 
is  the  total  thrust  produced  by  the  chaser  spacecraft  thrusters. 
Solving  for  the  control  bounds  for  the  defined  rendezvous 
scenario  and  the  assumed  chaser  parameters  gives  us:  - 
8.33  mm/s2  <  u  <  8.33  mm/s2 

The  above  described  optimal  control  problem  is  solved  nu¬ 
merically,  by  the  Gauss  Pseudospectral  Optimization  Solver 
(GPOPS)  [10].  The  GPOPS  is  an  open  source  code  that 
implements  a  direct  collocation  method  based  on  the  Gauss 
Pseudospectral  approach  for  solving  optimal  control  prob¬ 
lems. 

Propellant-optimal  trajectories  are  simulated  with  time- 
averaged  J2  with  corrections  model.  To  validate  the  correct¬ 
ness  of  the  derived  time-averaged  J2  equations,  the  trajec¬ 
tories  with  0°  inclination  are  simulated  and  compared  with 
Hill’s  equations  along  with  its  control  profiles.  It  should  be 
seen  that  for  this  specific  test  condition,  the  optimal  trajec¬ 
tories  in  both  case  should  be  perfectly  similar  along  with  its 
trust  profile.  Moreover,  the  effects  of  the  correction  errors  on 
the  optimal  trajectories  is  also  shown.  The  simulation  results 
obtained  including  the  optimal  control  profiles  are  discussed 
below. 


5.  Simulation  Results 

The  minimum-propellant-control  solution  with  the  optimal 
controls  defined  in  Equation  (24)  is  here  presented  first.  Ad¬ 
ditionally,  the  effect  of  the  optimal  trajectories  with  respect 
to  the  inclination  is  discussed,  it  also  shows  the  effect  of  J2 
at  different  inclinations  and  how  the  optimal  trajectories  are 
influenced  by  it.  Moreover,  Equation  (11)  which  incorporates 
only  the  differential  J2  without  corrections  is  studied  to  see 
how  the  change  in  the  reference  orbital  period  and  nodal 
drift  errors  affect  the  optimal  solution.  Finally,  optimization 
results  obtained  with  Hill’s  equations  only  are  presented  and 
compared. 

Minimum  Solution  for  Defined  Rendezvous  Scenario 

For  the  minimum-propellant-control  rendezvous  scenario  set 
in  Section  4  the  pseudospectral  method  yields  the  solution 
shown  in  Figure  [10].  Additional  planar  views  of  the  optimal 
trajectory  is  shown  is  figure  [11], 


Figure  10.  Optimal  trajectory  for  time-averaged  J2  with 
corrections  model  (iref  =  70°). 


The  green  marker  shows  the  initial  position  of  the  chaser  with 
respect  to  the  target.  The  final  maneuver  time  as  discussed  in 
the  previous  section  is  fixed  at  one  orbital  period  of  the  refer¬ 
ence  orbit  which  when  calculated  comes  to  tf  =  1.72  hr  and 


ends  at  the  red  marker.  It  is  to  be  noted  that  the  blue  markers 
on  the  trajectory  are  the  points  where  the  thrust  is  applied  and 
it  is  assumed  that  the  thrust  applied  is  impulsively. 


Figure  11.  Planar  view  of  the  optimal  trajectory  for  time- 
averaged  J2  with  corrections  model  (?Tef  =  70°). 


Figure  [12]  and  Figure  [13]  shows  the  time  history  of  the 
relative  position  and  velocity  of  the  optimal  trajectory  re¬ 
spectively.  Along  with  that,  the  control  profile  for  the  three 
normalized  components  of  the  thrust  force  is  displayed  in 
Figure  [14],  It  is  clearly  seen  that  as  there  is  no  significant 
cross-track  motion,  the  fz  control  is  not  very  active.  On  the 
other  hand  because  of  the  large  amplitude  of  the  motion  in  the 
radial  direction  (y)  and  the  two  rapid  consequential  change  in 
the  radial  directions  makes  the  fy  control  very  active  as  it 
maxes  out  on  the  thmst  bounds. 


Figure  12.  Position  profiles  for  time-averaged  J2  with 
corrections  model  (*ref  =  70°). 
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Figure  13.  Velocity  profiles  for  time-averaged  J2  with 
corrections  model  (zref  =  70°). 
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Figure  14.  Minimum  propellant  control  profiles  for  time- 
averaged  J2  with  corrections  model  (zref  =  70°). 


It  should  be  noted  that  the  relative  tolerance  of  10  3  was 
set  to  solve  the  differential  equations.  Addition  to  that  the 
terminal  tolerance  of  the  order  of  10-6  is  set  for  the  NLP 
solved  by  GPOPS,  which  does  not  necessarily  indicate  the 
quality  of  the  solution.  The  reason  is  that  for  the  solution  the 
parameters  of  the  trajectory  are  only  being  computed  at  41 
nodes  only  where  the  optimality  conditions  are  enforced. 

For  the  implementation  of  the  Gauss  Pseudospectral  method 
using  GPOPS-II  the  guess  for  the  initial  states  corresponds  to 
the  initial  conditions  and  the  guess  for  the  final  state  consisted 
of  zeros  for  the  entire  state.  The  guess  for  the  control  history 
was  zero  at  the  initial  and  final  times  for  all  controls.  The 
solution  obtained  appears  to  be  feasible  but  can  be  used 
only  for  offline  computations,  i.e.  on  an  open-loop  guidance 
scheme. 


Minimum  Propellant  Solution  for  Hill’s  Equation 

This  sections  draws  differences  in  the  optimality  between  the 
Hill’s  Equations  and  the  time-averaged  J2  with  corrections 
equation.  The  following  figure  shows  the  optimal  trajectory 
for  the  Hill’s  equations.  As  expected  a  solution  similar 
to  Hohmann  transfer  is  achieved.  It  can  be  hereby  seen 
that  the  use  of  the  Hill’s  equations  for  the  study  of  optimal 
trajectories  for  rendezvous  and  proximity  operations  presents 
very  conservative  results. 


Figure  15.  Optimal  trajectory  for  the  Hill’s  Equation  and 
time-averaged  J2  with  corrections  ( iref  =  0°). 


The  above  Figure  [15]  shows  the  propellant-optimal  trajec¬ 
tories  of  the  time-averaged  J2  with  correction  model  for  0° 
inclination  along  with  the  the  linearized  Hill’s  equations  of 
motion.  Comparing  the  plotted  optimal  trajectories  validates 
the  correctness  of  the  derived  time-averaged  J2  correction 
equations.  An  insignificant  error  noticed  in  the  relative 
position  is  in  the  same  order  of  magnitude  as  seen  in  Figure 
[8]  corresponding  to  the  curvilinear  coordinated  system  used 
for  the  derivation  of  the  time-averaged  J2  equations. 
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Figure  16.  Minimum  propellant  control  profiles  for  the  Hill’s 
Equation. 


Figure  [16]  and  Figure  [17]  shows  a  plot  of  the  resulting 
control  history  of  the  three  applied  thrust  forces  correspond¬ 
ing  to  the  optimal  trajectory  shown  in  Figure  [15]  using  the 
linearized  Hill’s  equations  of  motion  and  the  time-averaged 
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Figure  17.  Minimum  propellant  control  profiles  for  the  time- 
averaged  J2  with  corrections  model  (iref  =  0°). 


J2  with  corrections  model  for  0°  inclination.  Comparing  the 
identical  thrust  profiles  furthermore  validates  the  correctness 
of  the  derived  time-averaged  J2  equation  used  for  this  study. 

Figure  [18]  adds  to  the  study  by  showing  the  effect  of  the 
differential  J2  model  without  corrections  for  70°  inclina¬ 
tion  on  the  optimal  trajectory  profile.  This  simulated  orbit 
is  alongside  plotted  with  optimal  trajectories  of  the  time- 
averaged  J2  with  corrections  model  (iref  =  70°)  and  the 
Hill’s  equations  so  as  to  fully  understand  the  sizable  errors 
in  the  optimal  positions  caused  by  the  the  orbital  period 
mismatch  between  the  target  and  the  chaser  orbits  and  the 
nodal  drift  caused  by  the  separation  of  the  ascending  node. 


Figure  18.  Optimal  guidance  trajectory  for  different  models 

(iref  =  70°) 

Effect  of  Inclination  on  the  Perturbed  Optimal  Trajectory 

As  the  effect  of  J2  potential  is  directly  proportional  to  the 
inclination  of  the  target  of  the  reference  orbit.  A  generalized 
plot  below  in  Figure  [19]  that  shows  the  J2  perturbation  effect 
on  the  optimal  trajectories  as  the  inclination  of  the  reference 
orbit  is  increased. 

It  is  seen  that  at  lower  inclinations  the  optimal  route  of  the 
chaser  is  very  similar  to  that  of  the  Hill’s  equation  as  seen  in 
Figure  [15],  which  should  be  the  case  as  we  already  discussed 
that  the  effect  of  J2  is  not  very  significant  at  inclinations  near 
to  equator  (0°).  For  all  the  simulation  plotted  in  Figure  [19], 
the  same  rendezvous  scenario  is  considered  with  fixed  final 
maneuver  time  of  one  orbital  period. 


Figure  19.  Optimal  trajectory  for  time-averaged  J2  with 
corrections  model  for  different  inclination. 


6.  Conclusion 

In  this  paper,  the  optimal  close  range  rendezvous  problem 
for  a  ADR  spacecraft  to  approach  a  tumbling  upper  stage  for 
capture  is  formulated  and  analyzed.  A  set  of  new  linearized 
differential  equations  have  been  incorporated  for  describing 
the  relative  motion  of  satellites  in  the  presence  of  J2  potential. 
Effects  on  the  relative  position  of  the  satellite  due  to  time- 
averaged  J2  with  corrections  model  was  shown  to  advocate 
its  implementation  for  the  proposed  optimal  propellant  ren¬ 
dezvous  problem.  The  minimum-propellant-control  problem 
was  formulated  and  addressed  using  the  direct  collocation 
pseudospectral  methods.  Differential  J2  without  corrections 
model  was  all  analyzed  to  show  the  significant  errors  in 
the  relative  position  due  to  the  orbital  period  mismatch  and 
nodal  drifts.  The  desired  optimal  trajectory  of  the  chaser 
spacecraft  with  respect  to  the  tumbling  resident  space  object 
is  sought  such  that  the  desired  final  state  matches  its  position 
and  velocity  under  its  defined  tolerances.  The  effect  of 
J2  potential  is  examined  on  the  optimality  condition  for  a 
defined  rendezvous  scenario. 

Further  study  is  required  and  will  incorporate  the  the  full 
6dof  model  of  the  chaser  spacecraft,  where  both  the  rotational 
and  transnational  dynamics  will  be  considered.  Moreover,  to 
make  the  capture  scenario  more  realistic  collision-avoidance 
condition  will  be  imposed  in  form  of  a  path  constraint.  Along 
with  that  docking/capture-enabling  conditions  will  also  be 
incorporated  by  matching  the  chaser’s  and  target’s  docking- 
station  position  and  velocity  vectors.  Moreover,  optimal¬ 
time  and  optimal-energy  conditions  will  also  be  studies  for 
the  rendezvous  scenario  to  truly  find  the  optimal  trajectory 
for  target  capture.  A  nominal  validation  benchmark  will  be 
introduced  by  using  GMAT  (General  Mission  Analysis  Tool), 
an  open-source  pace  mission  design  tool  created  by  NASA 
Goddard  to  optimize  and  propagate  the  trajectories  under  the 
influence  of  J2  disturbance.  To  fully  incorporate  the  environ¬ 
mental  perturbations  that  affects  the  proximity  operations,  the 
authors  would  like  to  add  aerodynamic  drag  forces  into  the 
relative  equations  of  motion  model.  Its  combined  effect  along 
with  J2  geopotential  disturbance  will  be  studied.  Optimal 
trajectories  will  be  planned  with  all  the  above  constraints 
defined. 
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